Correlated singlet phase in the one-dimensional Hubbard-Holstein model 
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We show that a nearest-neighbor singlet phase results (Irom an effective Hamiltonian) for the 
one-dimensional Hubbard-Holstein model in the regime of strong electron-electron and electron- 
phonon interactions and under non-adiabatic conditions {t/too < 1). By mapping the system of 
nearest-neighbor singlets at a filling Np/N onto a hard-core-boson (HCB) t-V model at a filling 
Np/{N — Np), we demonstrate explicitly that superfiuidity and charge-density-wave (CDW) occur 
mutually exclusively with the diagonal long range order manifesting itself only at one-third filling. 
Furthermore, we also show that the Bose-Einstein condensate (BEC) occupation number no for the 
singlet phase, similar to the no for a HCB tight binding model, scales as -\/iV ; however, the coefficient 
of -\//V in the no for the interacting singlet phase is numerically demonstrated to be smaller. 

PACS numbers: 71.10.Fd, 74.20.-z, 71.45.Lr, 71.38.-k 
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I. INTRODUCTION 

The study of coexistence and competition between di- 
agonal long range orders [such as charge density wave 
(CDW) and spin density wave (SDW)] and off-diagonal 
long range orders (such as superfluid and supercon- 
ducting states) in electronic phases is a subject of im- 
mense ongoing focus. Of particular interest is the co- 
existence of CDW and superconductivity/superfluidity 
in layered dichalogenides (e.g., 2H-TaSe2, 2H-TaS2, and 
NbSe2)i, helium-4^, bismuthates (e.g., BaBiOs doped 
with K or P)^, quasi-one-dimensional trichalcogenide 
NbSej^ and doped spin ladder cuprate Sri4Cu2404i— , 
quarter-filled organic materials^, non-iron based pnic- 
tides (e.g., SrPt2As2)^, etc. 

Systems with more than one type of interaction typi- 
cally manifest a variety of phases of which some cooperate 
and some compete. A wealth of materials show evidence 
of strong electron-phonon (e-ph) interactions besides the 
ubiquitous electron-electron (e-e) interactions. For in- 
stance, transition metal oxides such as cuprates^iiS and 
manganiteaii"— and molecular solids such as fullerides^^ 
indicate strong e-ph coupling. The interplay of e-e and 
e-ph interactions in these correlated systems leads to co- 
existence of or competition between various phases such 
as superconductivity, CDW, SDW, etc. 

An archetypal model for understanding the co- 
occurring effects of e-e and e-ph interactions is the fol- 
lowing well known Hubbard-Holstein model (HHM)^^ 



where cj^ is the fermionic creation operator for itinerant 
spin-cr electrons with hopping integral t and number op- 
erator rij^ = c'j^Cja, Oj is the corresponding bosonic cre- 



ation operator characterized by a dispersionless phonon 
frequency cjq, with U and g representing the strengths of 
onsite e-e and e-ph interactions respectively. 

To understand the interplay between the e-e and 
e-ph interactions, the Hubbard-Holstein model has 
been extensively studied (in one-, two-, and infinite- 
dimensions and at various fillings) by employing vari- 
ous approaches such as exact diagonalizationi^"— , den- 
sity matrix renormalization group (DMRG )^^i^° , quan- 
tum Monte Carlo (QMC)^"— , semi-analytical slave bo- 
son approximations^i"— , dynamical mean field theory 
(DMFT)^"— , large-N expansion^, variational meth- 
ods based on Lang-Firsov transformatio n^^i^'^ , Gutzwiller 



approximatior 



and cluster approximatioi 



A6 



In our earlier worki^, in the regimes of strong Coulomb 
interaction and strong e-ph coupling, we derived an ef- 
fective Hamiltonian using a controlled analytic approach 
that takes into account dynamical quantum phonons. 
We solved this effective Hamiltonian numerically for fi- 
nite chains and presented a phase diagram for the one- 
dimensional Hubbard-Holstein model at quarter filling. 
It was shown in Ref. that while the e-e interac- 
tion produces nearest-neighbor (NN) spin antifcrromag- 
netic (AF) interactions which encourage singlet forma- 
tion, the e-ph interaction generates NN repulsion which 
is expected to promote CDW order. It was also shown 
that a correlated NN singlet phase occurs (at quarter- 
filling) and that it carries a signature of a CDW. In this 
paper, we demonstrate that the correlated singlet phase 
occurs at other fractions as well and analyze its nature. 
Our main result is the demonstration that the NN spin 
AF and NN repulsive interactions compete (instead of 
cooperate) to produce mutually exclusive (rather than 
coexisting) superfluid and CDW phases in the NN singlet 
phase. We show that the NN singlets manifest superfluid- 
ity (and no CDW) at all fillings that are less than one-half 
but not equal to one-third and a CDW state (and no su- 
perfluidity) at one-third flUing. Using a modified Lanczos 
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methodi^i^ and a newly developed world-line quantum 
Monte Carlo (WQMC) method we show that the singlet 
phase has no Bose-Einstein condensate (BEC) fraction. 

In the past, superconductivity due to onsite pairing 
has been a focus of a number of studies^"—. Here we are 
interested in a different situation, namely, NN singlets. 
Earlier a t-J-V model (involving bipolarons that are NN 
singlets) was introduced in Ref. [Sli This t-J-V model^ 
[that does not include the next-nearest-neighbor hopping 
terms but discusses them qualitatively] is similar to our 
effective Hamiltonian of Eq. (|9]) and can be regarded as 
a useful precedent and an endorsement of Eq. dH) . 

The paper is organised as follows: in Sec. H we briefly 
derive the effective Hamiltonian (that goes beyond the 
t — J model approximation of the Hubbard model by 
including the additional three site residue^"—) and ex- 
plain the various interaction terms and hopping terms. 
We also point out that the correlated singlet phase oc- 
curs at not only quarter-filling but also at other fillings. 
In Sec. HI, we show that the correlated singlet phase 
can be represented by a hard-core-boson (HCB) t—Vi—V2 
model. Next, in Sec. IV we discuss the possibility of 
formation of a CDW by mapping the t—Vi—V2 model 
onto the well understood t-V model. In Sec V, we obtain 
the superfluid density (in the thermodynamic limit) at 
different filling fractions by using finite size scaling. In 
Sec. VI, we analyze the BEC occupation number at vari- 
ous densities by employing the modified Lanczos method 
and a newly developed WQMC method. We close with 
concluding remarks in Sec. VII. 



II. EFFECTIVE HHM HAMILTONIAN 

We briefly outline below the procedure to get the 
effective Hubbard-Holstein Hamiltonian (with more de- 
tails being provided in Ref. [Tsh . Although we obtain the 
effective Hamiltonian here in one-dimension only, our ap- 
proach is easily extendable to higher dimensions as well. 
We first carry out the Lang-Firsov (LF) transformation^ 

^hh = e^Hhhe''^ where T = -5 ^jcriaj - a]) and 
get the following LF transformed Hamiltonian: 



JO- 



^ rij + {U - 2g^uJo) ^ n^^nj^, (2) 



where n^^ = d^^dj^- On dropping the last term, which 



is a constant polaronic energy, we recognize that Eq. Q 
essentially represents the Hubbard Model for composite 
fermions with Hubbard interaction Ueff = (U — 2g^uJo). 
In the limit of large Ueff/t, using standard treatment 
involving a canonical transformation, we get the following 
effective Hamiltonian written to second order in the small 
parameter t/Uefj^'—- 



t-J-ts 



;.d 



'j+laCLj_i„aja 



diir + H.c. 



JO- 



p. 



(4) 



4t^ 



, ^3 = J/4, S., is the 



where n'^ = n^^ + n^^, J = ^7=2^^ 
spin operator for a spin 1/2 fermion at site i, and Pg is 
the single-occupancy-subspace projection operator. Fur- 
thermore, the last two terms with coefficient t-^ (= J/4) 
arc the three site terms which when omitted from the 
above Hamiltonian Ht-j-t^ yield the well-known t — J 
Hamiltonian (for the composite fermionic operators dja)- 
The effective t — J — Hamiltonian, given in Eq. (U) , 
can be re- written in terms of the original fermionic oper- 
ators Cja as 



(5) 



where 

Hn = -te"»' 



Ci„ + H.c. P, 



^0 X! 



-J '^Ps {Sj ■ Sj+i 



p. 
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Ps 



Je-3 



Yl [4oCj+ioc]_iffCjff + H.c. 



Ps. (6) 



JO 



and 



where Xj = e^^""^ "j-* and Uj = rij^ -\- rij^. Next, we 
express as follows our LF transformed Hamiltonian in 
terms of the composite fermionic operator d^-^ = cj^^j: 

j 

+{U - 2g'uo) E n^^t^ji " ^'^o ^ (4t + "ji) ' (3) 
j j 



jo 



In the above equation, we have separated the Ht-j-ts 
Hamiltonian into (i) an electronic part Hq which is es- 
sentially a modified t — J ~ Hamiltonian containing a 
NN hopping with a reduced amplitude {te~^ ), electronic 
interaction terms with the same interaction strength J, 
three site terms with reduced amplitude Je~^ /4. and 
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no electron-phonon interaction; and (ii) the remaining 
perturbative part Hi which corresponds to the compos- 
ite fermion terms containing the e-ph interaction with 
= e^^^"^+i~°J^. Furthermore, since J/4 << t, we 
have ignored the following term in Hi 



Ps 



Je-9 



Je-3' 



J2 [c]^c,+iJ^_i^c,,{Z=^Zl - 1) +H.C. 



JO- 



[c],c,+iJ^_i^c,,{Z'^ZL - 1) + H. 



(8) 



where = e^3{a,-i-a^+i) ^ 

After carrying out perturbation theory to second-order 
(as outlined in Ref. and Appendix A), with ^/(gwo) 
as the small parameter—, we get the following effective 
Hamiltonian: 



H 



•ff 

hh 



'teffhti + Jhs — Vhnn — t2h„„ 

-{t2 + J^)Ka + J:iKs (9) 



where 



K ^Y.p^ (4+1.^^- + Ps, (10) 



hs ^^Ps (Sj ■ Sj+i - jn^Hj+i) Ps, (11 



hnn = ^(1 -nj+i5-)(l -njg){n.j„ - n^+icr)^, (12) 



JO 



and 



h„cr= ^(1 - nj + lg){l - njg){l - Uj-lg) 
JO 

X [c]-Hio(l ~ 2njCT)cj_iCT + H.c. , 
hcs= ^(1 - nj+ig){l - Jij-ia) 

JO 

X [c]oCj + loc]_i3Cjff -hH.c. , 

Ka= ^ ^lj + lo)(l - "jo)(l - nj^ig) 

JO 

X \^c]^Cj+iac]_i^Cjg +R.C. 



(13) 



(14) 



(15) 



The various coefficients are defined in terms of the sys- 
tem electron-phonon coupling g, the Hubbard interac- 
tion U, the hopping amplitude t, and the phonon fre- 
quency Wo as follows: V ~ /2g'^oJo, J ^ ^* 



teff = te-^\ t2 ~ t^e-sVs^wo, and J3 = Je-a' /i. 
Here the kinetic energy (which is small compared to 



the interaction energy) has contributions from four hop- 
ping terms: —teffhti corresponding to NN hopping (with 
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a reduced hopping integral teff= te~^ ), —t2h„„ repre- 
senting NNN hopping (with double-hopping coefficient 
t^e~^^ /q'^ujq), -{t2 + Jz)hag implying NN spin-pair 
aa hopping, and Jsh'^g leading to NN spin-pair aa hop- 
ping and flipping to aa; thus h'^g acting on a singlet 
state produces another singlet state, but with a nega- 
tive sign. The NN spin-spin interaction term Jhs (with 
-) and the NN repulsion term —Vhnn (with 



4t^ 



V ~ l^g^uj^) are the dominant terms in the effective 
Hamiltonian and compete to form a phase separated clus- 
ter at larger J (or smaller V jt at a fixed g and t/ujo). As 
J /V decreases, the cluster breaks up to undergo a dis- 
continuous transition to a correlated NN singlet phase as 
shown in the phase diagram [see Fig. [IJa)}^. At even 
lower values of J/V, we get separated single spins (rep- 
resented by isolated spin phase) with the transition at 
larger g being first-order while at smaller g it is weakly 
first order and not continuous [due to the fact that the 
system transforms from a superfluid to a CDW, i.e., tran- 
sition is between two phases of different symmetry}^. 
The prime objective of the current work is to character- 
ize the correlated singlet state. 

We will now compare the physics related to our ef- 
fective Hamiltonian, which accounts for various funda- 
mental processes involved in the kinetic and interaction 
terms, with the variational Lang-Firsov (LF) treatments 
reported^i^'^i^'^. As the degree of non-adiabaticity 
decreases, our NNN hopping term —t2haa contribution 
increases, effectively the hopping transport will be larger 

than that given by —teffht^^; these two hopping terms 

2 

together can be regarded as producing a less than e~^ 
suppression of the hopping integral reported in earlier 
variational LF treatments. Furthermore, concerning the 
effect of including a large Hubbard U term in a Hol- 
stein model; we get the NN interaction 2V reduced to 
2V — J/4; thus, the mobility would be enhanced which is 
consistent again with the earlier works using variational 
LF transformation. 



III. t-Vi-V2 HARD-CORE-BOSON (HCB) MODEL 

In the rest of the paper we study the correlated singlet 
phase. No pair of singlets can share a common site. The 
closest two singlets can approach each other is to have 
one spin from each singlet be on adjacent sites. The sin- 
glets transport via two processes: (i) the NN spin-pair aa 
hopping given by the h„g and h'^g terms in Eq. @ ; and 
(ii) a second order process involving breaking of a bound 
singlet state [with binding energy Eb = — J + i^/(g^a;o)] 
and hopping of the two constituent spins (of the singlet) 
to (a) neighboring sites in the same direction sequentially 
[yielding the term —t^h^g with tf, = t^e~^f /|i?_B|] or (b) 
neighboring sites in opposite direction and back [yielding 
the term —tbhnn]- Wc now make the important obser- 
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FIG. 1. (Color online) Plots obtained using modified Lanc- 
zos in a twelve-site system for t/uo = 1. Phase diagram in 
(a) depicts that the phase transition lines are close for both 
densities n = 1/4 and n = 1/6. Structure factor plots in (b) 
(drawn ai g — 2.2 and U/t = 17) for the effective Hubbard- 
Holstein model (HHM) of Eq. (O and the HCB t-Vi-V2 model 
of Eq. (|16|) showing that the two models are equivalent. 



vation that a NN singlet can be represented as a HCB 
located at the center of the singlet. Thus the system of 
NN singlets in a periodic lattice is transformed into a 
system of HCB also in a periodic lattice with the same 
lattice constant a but with the whole lattice displaced by 
a/2. Then the effective Hamiltonian of the HCB system 
is the following t-Vi-V2 model: 

j+i + H.c.) + ViUjUj+i + V'2»^in-j+2],(16) 

j 

where bj is the HCB destruction operator, rij = b'jbj, 
T = {t2 + 2J3 + tb), Vi = oo (because two singlets cannot 
share a site), and V2 ~ 2V - J/4 [with V2/T > 10 (i.e., 
V2/T >> 1) for parameter values in the singlet regime 
of our phase diagrams in Fig. [Ha)]. In the following 
we set T = 1. We corroborate our mapping of the ef- 
fective HHM Hamiltonian H^jJ (for the singlet phase) 
onto the HCB Hamiltonian Hi, by demonstrating in Fig. 
mb) that the static structure factor S{k) = J2i e*'''W^(0 
for the HHM and HCB cases coincide when the correla- 
tion function W{1) = {l/N)J2^[{A,A,+i) - {A,){A,+i)] 

is defined through Aj ee {S+ S^^^ + H.c.) for HHM and 
Aj = Uj for HCB. 

It should be made clear that, for performing calcula- 
tions, there is a distinct advantage of accessing bigger sys- 
tem sizes for the HCB system as compared to the HHM 
Hamiltonian. For instance calculations involving 8 HCB 
(equivalent to Sf and 84, electrons) on a 24 site lattice 
/ 24 \ 

require ( g ] = 735471 basis states in the occupation 

number representation and hence are certainly feasible 
using modified Lanczos method; on the other hand, using 
the same technique, one can barely deal with 8 electrons 
(4^ and 4^,) on a 16 site lattice for the HHM Hamiltonian 

as it requires I "'o^ 1 x ( " ) = 900900 basis states. It is 
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FIG. 2. (Color online) WQMC plot of the structure factor 
S{k) versus fc - for iV = L = 60, /3 = LAt with At = 
0.125, and at various densities - shows CDW at 71 = 1/3 
with S{Q) « A'^/9, i.e., maximum allowed value. The peak 
values S(Q) rapidly fall as n moves away from 1/3 and are 
independent of V2 at large values of V2 [see inset] . 
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FIG. 3. (Color online) Plots, obtained using WQMC at /3 = 
NAt with At = 0.125, showing correlations in the t-V\-V2 
model. The correlation function W{1), plotted for A'^ = 80 
sites in (a), does not seem to decay. The peak of the structure 
factor S{Q), plotted in (b) for various system sizes at n — 1/4, 
grows monotonically. 



also of interest to note that representing a NN singlet by 
a HCB located at the center of the singlet, although has 
been done here for a one-dimensional system, can also be 
done in higher dimensional systems. 



IV. CDW CORRELATIONS 

The repulsive terms in the HCB Hamiltonian Hh indi- 
cate that a CDW is possible. We study the correlations, 
by extending to our t-Vi-V2 model, the well documented 
WQMC approach for obtaining correlation functions and 
structure factor for the t-V model^. Plots of the struc- 
ture factor in Fig. [2] show a peak at wavevector Q = iirn 
suggesting a CDW. However (as shown in Fig. only at 
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filling n = 1/3, where the structure factor peak is approx- 
imately that for the strong CDW case corresponding to 
V2 ^ oo. can we assert that CDW occurs. Specifically at 
n ~ 1/3 and for V2 > 10, the W{1) has a simple structure 
[i.e., W{1) « 1/3-1/3x1/3 = 2/9 when Hs a multiple of 
3 whereas for other I values 1^(0 « -1/3 x 1/3 = -1/9] 
yielding S{k) « ^fc.2ir/3^/9- Furthermore (in Fig. [5]), 
the peak of the structure factor S{Q) (which remains es- 
sentially constant at all relevant interactions V2 > 10) 
rapidly decreases as n decreases from 1/3 - a trend that 
is similar to that of S{Q) for the t-V model as one moves 
away from half-filling^. Nevertheless, the plots of corre- 
lation function (in Fig. [3]) do not seem to decay at large 
distance (for both n = 1/4 and n = 1/5) while the struc- 
ture factor peak (for n — 1/4) seems to grow monoton- 
ically with system size - all indicative of a CDW. Later 
on, the above ambivalence will be resolved and it will be 
demonstrated unequivocally that our t-Vi-V2 model has 
a CDW only at n = 1/3 while at other fillings n < 1/3 
superfluidity (and no CDW) results. 

Since Vi = 00 and because we are dealing with a 
one-dimensional system, we simplify the phase transition 
analysis by performing an exact mapping of the TV-site t- 
V1-V2 model onto a t-V model with N — Np sites and with 
V ^¥2- This enables us to access bigger system sizes for 
performing numerics; furthermore, since the phase di- 
agram of the t-V model is well known, we can clearly 
determine the existence of a CDW which was not possi- 
ble using the above structure- factor/correlation- function 
analysis. Later, we will also show that the t-V model 
lends itself to a simple finite size scaling approach for 
obtaining accurately the superfluid density in the ther- 
modynamic limit. 

We first recognize that we can recast the HCB Hamil- 
tonian in Eq. (|16|) as the following projected Hamiltonian 
where NN sites of a particle are projected out: 

Hb =E[-^i(l - n,-i)6]5,+i(l - n,+2) + H.c.} 
j 

+ V2{1 - nj_i)?ij(l - nj+i)7ij+2(l - nj+3)] 
=E[-T(6]6,+i + H.c.) + V2n,fi,+2]. (17) 



0.95 



where 6^ = (1 — nj_i)6j(l — 7ij+i) and fij = b^jbj. Next, 
we observe that commutes with 'Y^- (1 — nj+i) and 
thus the total number of excitons (with each exciton com- 
prising of a particle with a hole to its right) is conserved. 
Physically, this is due to the fact that infinite NN repul- 
sion ensures that the neighboring sites of a particle are 
unoccupied. With each particle, we associate only one 
neighboring vacant site (say, the site on the right side of 
the particle) so that situations such as particles on NNN 
sites can also be dealt with. Then by deleting the sites 
of the holes in all the excitons and having only a NN in- 
teraction V = V2 and no other interaction in the reduced 
system of Ni = N — Np sites, we get the same eigenener- 
gies (see Rcf. [gO for a similar analysis for the t-V model 
in one-dimension). We further recognize that there is a 
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FIG. 4. (Color online) Superfluid density for an infinite sys- 
tem n*'' at various densities n and interactions V2 for the 
t-Vi-V2 model at Vi = 00 are depicted in (a). Values of n^'' 
in (a) are the intercepts, obtained by extrapolation of the 
straight lines through the tIs data plotted at various 1/Ni 
values, in figures such as (b) and (c). The solid lines in (a) 
are for V2 = 00 and obtained from Eq. (f20)) . 



one-to-one mapping between the eigenstates of the H^^ 
Hamiltonian and the eigenstates of the t-V Hamiltonian 
Ht~v, 

Ht-v = E[-T(6]&j+i + H.c.) + Vn,n,+i], (18) 



with V = V2 and Ni sites while the corresponding 
eigenenergies are identical. We can thus extract the 
eigenenergy spectrum of the t-Vi-V2 model by study- 
ing the equivalent t-V model. We first observe that 
n = Np/N ~ 1/3 for the t-Vi-V2 model corresponds to 
the n = Np/{N - Np) = 1/2 for the t-V model and 
thus superfluid density vanishes (as the two models have 
the same eigenenergies) and a CDW results^ since the 
mass gap is the same for both. Furthermore, at all frac- 
tions 71 < 1/3 for the t-Vi-V2 model we get a superfluid 
(and no CDW) since for the t-V model the same is true 
at n < 1/3^. Lastly, since n = 1 for the t-V model 
translates to n = 1/2 for the t-Vi-V2 model, we note 
that electron-hole symmetry for the t-V model guaran- 
tees that t-Vi-V2 model exhibits superfluidity and ab- 
sence of CDW for 1/3 < n < 1/2 as well. 



V. SUPERFLUID DENSITY 

We will now substantiate the above observations on the 
occurrence of superfluidity through calculating the super- 
fiuid density by threading the chain with an infinitesimal 
magnetic flux. We will exploit the one-dimensionality of 
the system and outline a simple flnite size scaling ap- 
proach to calculate the superfluid density in the thermo- 
dynamic limit. We first note that the energy for the t- 
V1-V2 model, when V2 = 00 and (as before) Vi = 00, 
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FIG. 5. (Color online) Superfluid density decaying exponen- 
tially with system size for the CDW state at one-third filling 
and large NNN repulsion V2. 



is given by the tight binding Hamiltonian energy for 
N2 = N — 2Np particles where we have excluded both 
the NN and NNN holes to the right of the particles in 
the t-Vi-V2 model. The total energy, when threaded by 
a flux 0, is expressed as 



E{9) = -2r^cos[fc + e/N2 



(19) 
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FIG. 6. (Color online) Plots of BEC occupation number no, 
obtained from modified Lanczos (open circles) and WQMC 
(crosses), with (a), (c), and (e) pertaining to t-Vi-V2 model 
(with Vi = CO, and V2 = 35) while (b), (d), and (f) respec- 
tively pertaining to the corresponding tight binding model 
with enhanced densities Np/{N - 2Np). For WQMC, /3 = 
NAt with At = 0.125, 0.15, and 0.175 for (a), (c), and (e) 
respectively. 



Then the superfluid fraction is given by*^' 



6132 



NpT 



2 9612 



sin 



ttNj_ 

Nn 



9=0 



N, 



p sm 



(20) 



where anti-periodic (periodic) boundary conditions have 
been taken for even (odd) values of Np. The superfluid 
density in the thermodynamic limit n*,^ can be related 
to the finite (A'^2-site) system superfluid density rij as 
follows: 



Ah 









f ^ \ 


I 


KN2) ^ 


120 


{N2) 



(21) 



From the above expression (valid for V2 = 00), at a flxed 
density, we expect (n*'* — ns)/ns (x l/iV| or l/Nl (with 
corrections of order 1 /N2 or 1/Nf) for the large but finite 
V2 case as well. We calculated the superfluid density at 
various large values of V2, system sizes N, and fllling 
fractions n; we find [as exemplified in Figs. Sljb) andSlJe)] 
that Us indeed varies linearly with 1 /N^ using which we 
obtain the various rt*'* values. 

From Fig. Hl^a), we see that the superfluid density 
(plotted in the thermodynamic limit) gradually decreases 
with increasing V2 and reaches the asymptotic value; 
the nl^ values for smaller filling fractions decrease more 
slowly because repulsion is less efl'ective at lower den- 
sities. Regarding the superfluid density at n = 1/3 and 
V2 = 00, it vanishes at all system sizes as can be seen from 
Eq. ((20)) . However, at finite V2 > 10, Ug vanishes expo- 
nentially with system size [as shown in Fig. which is 
consistent with the fact that there is a full CDW gap at 
n = 1/3. 



VI. BEC OCCUPATION NUMBER 

Lastly, we will calculate the Bose-Einstein condensate 
(BEC) occupation number uq. We first recall the well- 
established result that no, for a system of HCB in a one- 
dimensional tight binding lattice, varies as C{n)^/N in 
the thermodynamic limit with the coefficient C (n) mono- 
tonically increasing from as the density n increases 
from to 1/3^1^; consequently, the condensate fraction 
no/Np cx 1/\/N — 7> 0. Next, in the presence of repul- 
sion (as argued below), we expect the BEC occupation 
number tiq to again scale as \/N ; however, the coefiicient 
of y/N will be smaller due to the restriction on hopping 
imposed by repulsion. 

The Bose-Einstein condensate (BEC) occupation num- 
ber no is obtained from 



(22) 



where I'i'o) is the ground state. We calculate uq using 
two methods - modified Lanczos for smaller systems and 
a newly developed WQMC method for both small and 
larger systems (see Fig. |6]). The values of no for our t-Vi- 
V2 model in a iV-sitc original system So at various den- 
sities [such as n = 1/4, 1/5, 1/6] seem to be smaller than 
the no for the corresponding transformed tight binding 
system S2Np, realized when Vi = V2 = 00, with N — 2Np 
sites and enhanced densities [n/(l — 2n) = 1/2, 1/3, 1/4, 
respectively]. This can be understood from the fact that, 
in the transformed S2Np system of iV — 2Np sites [based 
on Eq. (HH)], a particle can hop to more sites between 
two particles than in the original t-Vi-V2 system leading 
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to a larger uq. For the S2Np system, it is important to 
realize that rto oc y^N — 2Np cx -y/iV- 

We will now consider a tight binding system S^Np with 
N — ANp sites and Np particles so as to obtain the lower 
bound for the BEC occupation number hq for the A^- 
site t ~ Vi ^ V2 system So- For every configuration in 
the S^N^ system, there is a corresponding configuration 
in the So system that can be obtained by adding two 
empty sites to the right and two empty sites to the left 
of all particles. Furthermore, the ground state kinetic 
energy contribution of the SiNp and S2Np systems are 
both proportional to N; hence, in the ground state of the 
original So system, the combined probability weighting 
of all the configurations obtained from the S^Np system 
(by adding 4 empty sites next to every particle) is a finite 
fraction. Since the BEC occupation number hq of S^Np 
system scales as a/TV, it follows that the lower bound 
of the hq for the original So system also varies as \/N. 
Thus, the BEC occupation number ng of the original 
A^-site t—Vi — V2 system So will vary as \/N since it is 
constrained from above by uq oc -v/ZV for the 5*2^^ system. 

At higher densities (i.e., 1/3 > n > 1/5) in our t-Vi- 
V2 model, we find that the values of no seem to increase 
more slowly with system size [see Figs. in{a),[nic), and 
[S]Je)] - this being due to smaller coefficients of ^/N re- 
sulting from interaction effects. Moreover, we also note 
[from Figs, ^h) and|6l[e)] that the value of no [i.e., the 
coefficient of ^/N in the expression for no] decreases due 
to repulsion. 

Our new WQMC method (see Appendix B for details) 
to obtain BEC fraction is a modification of the standard 
approach to studying correlations in the xxz model.— 
Since the Hamiltonian is real, it can be shown that the 
probability amplitude of any basis state in the ground 
state expression can be taken as real and non-negative. 
Consequently, we approximate the ground state by 

|*)=^yiSeME|^^,), (23) 

i 

with Z being the partition function, \(j>i) a basis state 
of the system in the occupation number representation, 
and (3 being sufficiently large. Then we calculate no by 
setting l^-o) = [*) in Eq. (HH). Our WQMC approach 
to no has been benchmarked against the modified 
Lanczos method for small system sizes (see Fig. 
The number of passes needed to estimate l^*) turns out 
to be an order of magnitude larger than that needed 
for obtaining correlation functions by WQMC. We 
take I^P) to be the state that produces an estimate of 
the kinetic energy (^'|A'|^') (with K being the kinetic 
energy operator) that is closest to the usual WQMC 
estimate \ exp[- pH]K\(j>i) / \ exp[- I3H]\(I)i))qmc 

where ()qmc denotes a quantum Monte Carlo average 
over various states \(t>i). 



VII. CONCLUSIONS 

In this paper, we have analyzed the correlated NN 
singlet phase predicted by the effective Hamiltonian of 
the Hubbard-Holstein model by essentially mapping the 
Hamiltonian onto the well-understood one-dimensional 
t-V model with large repulsion. Because the physics is 
dictated by the t-V model, we find that CDW and super- 
fluidity occur mutually exclusively with CDW resulting 
only at n = 1/3 while superfluidity manifests itself at all 
other fillings. We also show that the the BEC occupa- 
tion number no for our model scales as \/N similar to 
the no for a HCB tight binding model; additionally, we 
demonstrate numerically (using a new WQMC method 
and a modified Lanczos algorithm), at n ^ 1/3, that the 
no for our model is smaller than the no for a HCB tight 
binding model. 

We close by observing that, while CDW and supercon- 
ductivity seem to be incompatible in the one-dimensional 
HHM, experimental results (such as those reported in 
Refs. [ll-Q) suggest that they can coexist in higher dimen- 
sions. Furthermore, the vanishing of BEC fraction for the 
HHM is again an artifact of the one-dimensionality and 
should make way to non-zero fractions for higher dimen- 
sions just as in the case of the xxz model^. 
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Appendix A 

In this appendix, we will outline our approach to 
carrying out perturbation theory and obtaining the 
ground state energy. We assume a Hamiltonian of the 
form H = Ho -\- Hi where the unperturbed Ho has 
separable eigenstates \n,m) = \n)ei ® \m)ph with |0, 0) 
being the ground state with zero phonons; the eigenen- 
ergies, corresponding to |n, m), are En}n = E'^ + E^J^. 
Furthermore, the perturbation Hi is the elcctron-phonon 
interaction term of the form given in Eq. ([7]) . 

After a canonical transformatioi*i^, we obtain 

H = e^He-^ 

=Ho + Hi + [Ho + Hi,S]+^[[Ho + Hi,SlS].iAl) 

In the ground state energy, we know that the first- 
order perturbation term is zero by construction (in fact, 
(ni, 0|i?i|n2, 0) = 0). To eliminate the first-order term 
in Hi, we set Hi + [Hq, S] ~ 0. Consequently, we obtain 
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the matrix elements 



(ni,mi|5'|n2,m2) = - 



(ni,mi|iJi|n2,m2) 



(A2) 



We now assume that both NN hopping integral te~^ 



J 



and the Heisenberg spin interaction strength J are much 
smaller compared to the phononic energy ujq which is true 
at large couplings g. Hence, we make the approximation 
{Elvira,- E^^l^,) ~ {EP^^-EP^J; then, using Eqs. dS!]) 
and (IA2I). we obtain 



ph{mi\H\m2)ph — ph{mi\Ho\m2)ph + t:'^ ph{mi\Hi\m)ph ph{m\Hi\m2)ph 



pph pph rpph j-^ph 
rjm2 -t^m -C^mi -C/m 



(A3) 



(2) 

Next, it is important to note that the second order correction En,m, corresponding to the unperturbed eigenenergy 
En}ni, can be expressed as follows: 



(2) 



E 



{n,m\Hi\m)ph ph{m\Hi\n,m) 



(A4) 



Furthermore, since (rii, 0|-ffi|n2, 0) — 0, (n, 0|iJ|n, 0) is 
the total energy that resulted from performing second or- 
der perturbation theory on the unperturbed energy e'^^'^q. 
Our procedure for finding ground state amounts to ob- 
taining the lowest eigenvalue for the matrix with elements 
(m, 0|if|n2, 0); this is equivalent to finding the ground 
state of the effective Hamiltonian TJg (as was done in 
Ref. [H): 



where 



p/i(0|i7i|m)p/i X p/i(m|ffi|0)p/,. 



Et - E^ 



(A5) 



(A6) 



This procedure amounts to considering the restricted 
subspace spanned by eigenstates |n, 0)i obtained from 
carrying out first order perturbation theory on |n, 0): 



,0)l = |71,0)+^ 



\m)ph ph{m\Hi\n,0) 



ph 



(A7) 



It is important to recognize that the state In, 0)i is not 
separable, i.e., cannot be expressed as a product of an 
electronic wavefunction and a phononic wavefunction. 
We have restricted ourselves to the subspace of the states 
|n,0)i because the states \n,m ^ 0)i correspond to 
higher energy states due to the fact that the electronic 
excitation energy is much smaller than the phononic en- 
ergy, i.e., te~^ « ujq. Additionally, we would like to 
point out that the total ground state energy (in second 
order perturbation theory) is obtained by diagonalizing 
the matrix whose elements are (ni, 0|iJ|n2, 0)i. 



Appendix B: WQMC FOR BEC FRACTION 

We will discuss, in brief, the usual world-line quan- 
tum Monte Carlo (WQMC) approacb^i^ adapted for 



r 



calculating correlations in our t-Vi-V2 model Hamilto- 
nian given below: 

H, =Y.H^= Y}-nh%+i + H.c.) 
J J 

+ VirijUj+i + V2rijrij+2]. (Bl) 

Since this is quite similar to the t-V model, we can em- 
ploy the checkerboard decomposition Hb = Hi + H2 
where iJi = ^ and H2 = J2 ■ 1^ impor- 

j odd j even 

tant to note that both Hi and H2 consist of independent 
two-site pieces. Because of the decomposition, it becomes 
easier to evaluate the expectation value of an operator A 
given by 



{A) 



(B2) 



with A involving only number operators (such as n.inj) 
or NN hopping operators (such as 6]6_,+i + H.c). Now 
we calculate the partition function: 

Z = Trie-^""] 

i2L\U2\i2L-l) ■■■{h\Ul\i2) {i2\U2\il) ■ 

ix.....i-2L 

Here Ui ~ e^^'^^\f3 = LAr, and each of \ii), ...,\i2L) 
form a complete basis set in the occupation number rep- 
resentation. Here the world lines are the locus of the 
particles in the imaginary time (r) direction. 

For the density-density correlation function {niUi^i) 
(which is the expectation value of a diagonal operator), 
the above procedure of inserting 2L time slices yields the 
simple form 



1 



{iL\nini+i\iL) + {iL+i\nini+i\iL+i)]) 



QMC ■ 



where ( )qmc represents average over many QMC 
passes. Notice that we have concentrated only on L and 
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L + I time slice indexes although expectation value can 
be taken over all the 2L time slice indexes for better 
statistics. As for (6jfej+i +H.c.) (which corresponds to a 
non-diagonal operator), WQMC procedure yields 



{b]bj+i + H.c.) = 



{iM\{blbj+i +R.c.)Uk\iM+i 
(u/|C/fc|iM+i) 



)QMC, 



where, for odd (even) values of j, we take fc = 1 (2) 
and even (odd) M. However, as regards obtaining ex- 
pectation value of {bjbj+m + H.c.) for m > 1, the simple 
procedure (involving checkerboard decomposition) given 
above is not applicable; moreover, other suggested pro- 
cedures in the literature are complicated^. 

Here, we propose an alternate simple method for eval- 
uating (6j6j-f,„ -|- H.c.) for TO > 1 and thus obtaining the 
BEC occupation number 



(B3) 



with jvpo) being the ground state. To the WQMC method 
mentioned above, we add our trick to construct jv^o) as 
a linear combination of the basis states |(/)i) in the occu- 
pation number representation, i.e., |\E'o) = '^ailcpi) with 



1. Once we get a good estimate of the ground 



state l^'o), we can calculate the expectation values of any 
operator. 

After equilibrium (which is attained after several QMC 
passes), we run the simulation for a sufficient number of 
QMC passes and store the basis states corresponding to 
time slices L and L -I- 1 in each pass. It is obvious that 
some of the basis states will occur more frequently. The 
frequency of occurrence of a basis state \(f>i) is propor- 
tional to the probability (of) of its occurrence in the ex- 
pansion of the ground state jvPo)- Now, the coefficients 
Oi can be taken as real because the Hamiltonian is real 
and consequently I'J'o) can also be taken as real. Further- 
more, all a,; can be taken to be positive for the following 



reason. Firstly, the expectation values of NN and NNN 
interaction terms remain unaffected by the sign of Oi . 
Next, the expectation value of the hopping term is given 
by 

-T{^o\blbi+M = -T^(0,|a.(&|6i+i)a,|</),)] 
= -T'^{(j)i\aiCk\4'k)] 

i,k 

= -TJ2a^c^. (B4) 



This value is minimized when and Ci have the same 
sign. Then, if we take a.; to be positive for all i, > 
for all i. Thus in |5'o) — J^^il't'i): we can take all to 

i 

be positive and real. 

Let and Ej be the eigenstates and the eigenener- 
gics of the Hamiltonian with Eq being the ground state 
energy. For sufficiently large /3, we approximate the 
ground state by 



(0,|exp[-/3iJ]|(/., 



(B5) 



because then 



i 




E,|*,)(*,|exp[-/?F]Efcl*fc)(* 



(0,|^'o)cxp[-/?£;o](*o|0,) 



z 



\<f>i 



(B6) 



since the partition function Z = J2i{'^i\ exp[— 
exp[-~l3Eo]. 
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